EFFECTIVE METHOD OF COMPUTING LI'S 
COEFFICIENTS AND THEIR PROPERTIES 



KRZYSZTOF MASLANKA 

Abstract. In this paper we present an effective method for computing 
certain real coefficients A„ which appear in a criterion for the Riemann 
hypothesis proved by Xian-Jin Li. With the use of this method a se- 
quence of over three-thousand A n 's has been calculated. This sequence 
reveals a peculiar and unexpected behavior: it can be split into a strictly 
growing trend and some tiny oscillations superimposed on this trend. 



1. Introduction 

Since its formulation almost a century and a half ago the Riemann hy- 
pothesis (hereafter called RH) is commonly regarded as both the most chal- 
lenging and the most difficult task in number theory [14] . It states that all 
complex zeroes of the zeta function, defined by the following series if 5Rs > 1 



CO 



n° 

n=l 



and by analytic continuation to the whole plane, are located right on the 
critical line = |. RH, if true, would shed more light on our knowledge 
of the distribution of prime numbers. More precisely, the absence of zeroes 
of £(s) in the half-plane $ts > 8 implies that (see JHj, theorem 30) 

(1.2) it(x) = li(x) + 0(x d logx) 

where ir(x) is the number of primes not exceeding x and li(x) denotes loga- 
rithmic integral. Therefore, the value 6 = \ ( as Riemann conjectured) makes 
the theorem useful since the error term in p.2[) is the smallest possible. We 
do know that on the critical line lie infinitely many complex zeroes [3] and 
that among several billions of initial zeroes there is no counterexample to 
RH, cf. [EH, Q7]. 



(Place Figure 1 about here.) 



Date: 9 May 2004. 

Key words and phrases. Riemann zeta function, Riemann hypothesis, Li's criterion, 
numerical methods in analytic number theory, 
e-mail: maslanka@oa.uj.edu.pl. 

1 



2 



KRZYSZTOF MASLANKA 



2. Li's Criterion 



In 1997 Xian-Jin Li |1U| presented an interesting criterion equivalent to 
the Riemann hypothesis: 

Theorem 2.1. RH is true if and only if all coefficients 

1 dP 



(2.1) 



A n 



8=1 



T(n) ds r 

are non-negative, where 
(2.2) 

An equivalent definition of A n is (see jTU], formula 1.4): 
(2.3) A n = Y, ( 1 



e( s ) = 2( s -i)7r-/ 2 r (i + £)c< 



1-i 



where the sum runs over all (paired) complex zeroes of the Riemann zeta- 
function. However, the above definitions of A n are not suitable for numerical 
calculations. In this paper I shall present an effective method for calculating 
these coefficients. The gathered data investigated numerically up to n = 
3300 reveals unexpected properties: it contains a strictly growing trend plus 
extremely small oscillations superimposed on this trend. 

The following decomposition of A ra is implicitly given in a recent paper by 
Bombieri and Lagarias ([2], Theorem 2): 



(2.4) A n 



l-(log(4vr)+ 7 )- + ^(-l) J 
i=2 



(1-2-*) CO") 



trend 



-E 

3=1 



n 



oscillations 
= A n + A n 

Using the language of signal theory (perhaps not very common but some- 
times appropriate in number theory) one can say that the decomposition 
(|2.4j) uniquely "splits" the behavior of the sequence of {A n } into a strictly 

growing trend A n and certain tiny oscillations X n superimposed on it. It 
may be proved that the trend can be expressed as 

1 d n 



(2.5) 

while 
(2.6) 



At. 



r (n) ds r 
while the oscillations are 



„n-l 



in(^/ 2 rfi + ^ 



A„ 



1 



T (n) ds 



n 

-=-[ a -lln(( a -lK«)] 



s=l 
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It may also be proved that the trend is indeed strictly growing as n tends 
to infinity. It is evident that 1)2.5)) differs from the main definition 1)2.1)1 
simply by replacing £(s) by the much simpler function 7r _S//2 r(l + s/2). On 

the other hand, the oscillatory behavior of A n is not so evident, nevertheless 
it may be investigated numerically. I shall return to this decomposition 
later. 

The problem of calculating numerically both components of \ n (i.e. trend 
and oscillations), using directly 1)2.5)1 and 1)2.6)1 . is rather hopeless, not to say: 
malicious. One has to take the n th derivatives of functions which depend of 
variable s and are labelled by parameter re. When n tends to infinity both 
families of differentiated functions tend to right angle shaped figures and the 
derivatives are to be taken just at the almost singular point s = 1. What 
is interesting, the first derivative in s = 1 for all functions related to the 
oscillating part is the same and equal to the Euler constant. 

However, it is possible to calculate several tens of initial derivatives using 
direct numerical approach, for example Mathametica's ND built-in function. 
This function has several parameters which enable to control the required 
accuracy. One must be aware, however, that there is no guarantee that the 
result will be correct. As with all numerical techniques for evaluating the 
infinite via finite samplings, it is sometimes possible to "fool" ND into giving 
an incorrect result. 

(Place Figure 2 about here.) 

It turns out that the following function fits very well to the numerically 
tabulated values of 1)2.5)1 : 

(2.7) a (1 + rehire) + cn 

with 

a = - ± 8 ■ 1(T 9 
2 

c = -1.130330701... 

(A choice very similar to 1)2.7)1 was suggested to me by J. Lagarias, [S].) 
Recently I learned that A. Voros JH] using classic technique of saddle-point 
method calculated the exact value of c 

c = — (7 — 1 — In 2"7r) 

where 7 is Euler constant. Simple fitting procedure gives also coefficients of 
consecutive terms which appear to be related to Bernoulli numbers B k 

B k 1 1 t 1 1 | 1 

~ 2^ 2kn k ~ 1 ~ 4 ~ 24ra~ + 240n 3 ~ 504n 5 + 480n 7 

k=l 

This series works well for a dozen or so initial terms although it is formally 
divergent. 
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The starting point of Li's approach to RH is a certain transformation of 
the complex plane into itself using the map s h z = 1 — 1/s (which is a 
special case of Mobius transformation). Under this transformation the half- 
plane > g is mapped into the unit disk (with the critical line = | 
becoming the unit circle, see Figs. 3 and 4). This was Li's original idea. 
However, he was inspired by studying A. Weil's proof of RH for function 
fields over finite fields where the critical line is transformed into a unit circle 

nu. 

(Place Figures 3 and 4 about here.) 



3. The main derivation 

It has been known since medieval times that the harmonic series 1 1/k 
diverges. This was proved long ago with the use of elementary methods by 
Nicole d'Oresme in the 14th century, and, much later, independently by 
Pietro Mengoli (in his book on arithmetic series Novae quadraturae arith- 
meticae, 1650) as well as, using yet another method, by the Bernoulli broth- 
ers. 

A natural question emerges: how fast does this series diverge? It turns out 
that its divergence is "weak" , more precisely: logarithmic. The quantitative 
answer to this question implies the definition of the following famous number 
called the Euler-Mascheroni constant: 



:= lim f - 

\k<x 



(3.1) 7 := lim | > - - logx | = 0, 5772156649... 



Its natural generalization is the sequence j n defined by 

vn (log*)" +1 



( _-\\n I 1 

(3.2) ln-= K -J- lim £r(logfc)' 

\k<x 



n+1 



where 7 = 7- These are the so-called Stieltjes constants 1 . Another "simi- 
lar" very useful sequence denoted by rj n is defined by 



(3.3) rj n :-- 



n\ x^oo \ k n + 1 

\k<x 



(logx) 



^It should be noted that the function StieltjesGamma[n] implemented in Wolfram's 
Mathematica, which employs Keiper's algorithm |7j, uses a different convention. It is re- 
lated to our 7 n via 

(-1)™ 

7 n = - — ■— StieltjesGamma[n] 
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where A(fc) is the so-called von Mangoldt function defined for any positive 
integer k as: 



(3.4) A(k) 



log p if k is a prime p or any power of a prime p n 
otherwise 



The above sequences are important on their own right since they appear 
in the Laurent expansions for £(s) and its logarithmic derivative around 
s = 1. (There are different conventions when defining these numbers, here I 
have adopted those of Bombieri and Lagarias [2] ) : 



1 oo 

(3.5) C(- + l) = - + E 



7n« 
ra=0 



oo 



(3.6) J' {s + 1) = ± + ^ VnS n 

* n=0 

Integrating the second equation (|3,6|) with respect to s, inserting the result 
into the first one and equating coefficients in the appropriate powers of the 
variable s one can find explicit relations between the 7 n and the rj n : 



n+1 



( 3 - 7 ) 5>NTXT = -!°g 1 + 

n=0 \ n=0 / 

00 Q .n+1 00 ( i \k / oo \ k 

E%^T - E^'E^" 

n=0 fc=l \n=0 / 

(k) 

Now introduce the coefficients ch defined by 

oo / oo \ k 

n=0 \n=0 / 

Employing a certain formula from [1] (formula 0.314, i.e., raising a power 
series to an arbitrary integral exponent) one can express the c coefficients 
by the following recurrence relations: 



(3.8) 



"0 — 7 

1 m— 1 

$ = — E[ fcm -( fc+1 )^--^ fe) 

7717 — ^ 
' i=0 
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The matrix of coefficients c depends on {7 n }: 

m = m = 1 m = 2 m = 3 m = 4 



/c = 1 



7o 7l 72 73 74 

„2 o„, „, „,2 



(3.9) k = 2 7 g 2 7o7l 7 f + 2 7o72 2 7l72 + 2 7o73 
k = 3 7 jj 3 7 g 7l 3 7o7 f + 3 7 g 72 
A; = 4 7 <* 4 7 jj 7l 
A; = 5 t§ 

(In what follows only the upper triangular part of this infinite matrix will 
be needed.) 

(Place Figure 5 about here.) 



With the help of (|3.7[) the coefficients r\ n may further be expressed using 
the elements of the matrix c as 

\k+l 



(fc+i) 

k 

k=0 



From this we have: 
(3.11) rj = - 7o 

Vi = +7o- 2 7i 

V2 = -7o + 3 7o7i - 3 72 

V3 = +7o - 4 7o7i + 2 7 ? + 4 7o72 - 4 7s 

ri4 = -7o + 5 7o7i - 5 7o 7? - 5 7 q7 2 + 5 7i7 2 + 5 7 7 3 ~ 5 7 4 > 



Finally, the oscillating parts of A ra are expressible as polynomials in the 
Stieltjes constants: 

(3-12) ^ = -E("j^-i- 

Using now (|3.11|) and (|3.12|) we finally obtain: 

(3.13) Ax = 7o 

A2 = 2 7o - 7 q + 2 7l 

A3 = 3 7o - 3 7 q + 7 q + 6 7l - 3 7o7l + 3 72 

A 4 = 4 7o - fryg + 4jI - ^ + 12 7l - 12 7o7l + 4 7 § 7l - 2 7 ? 
+12 72 - 4 7o72 + 4 7s 
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Here are some numerical values of various numbers used in this paper: 



n 


7n 


Vn 


\ 


An 





+0.577215664902 


-0.577215664902 


' 


" 


1 


-0.0728158454837 


+0.187546232840 


0.577215664902 


0.0230957089661 


2 


-0.00969036319287 


-0.0516886320332 


0.966885096963 


0.0923457352280 


3 


+0.00205383442030 


+0.0147516588255 


1.22069692822 


0.207638920554 


4 


+0.00232537006547 


-0.00452447788850 


1.37558813187 


0.368790479492 


5 


+0.000793323817301 


+0.00144679520453 


1.45826850020 


0.575542714461 


6 


-0.000238769345430 


-0.000471544078185 


1.48829832721 


0.827566012282 


7 


-0.000527289567058 


+0.000155180294164 


1.48019084024 


1.12446011757 


8 


-0.000352123353803 


-0.0000513452121181 


1.44485574412 


1.46575567715 


9 


-0.0000343947744181 


+0.0000170413570471 


1.39059640679 


1.85091604838 


10 


+0.000205332814909 


-5.66605092104-10 -6 


1.32380368370 


2.27933936319 


11 


+0.000270184439544 


+ 1.88584861186-10~ 6 


1.24944277582 


2.75036083822 


12 


+0.000167272912105 


-6.28055422786 10~ 7 


1.17139824694 


3.26325532062 


13 


-0.0000274638066038 


+2.09240519074-10~ 7 


1.09272131711 


3.81724005785 


14 


-0.000209209262059 


-6.97247031237 10 -8 


1.01580941259 


4.41147767868 


15 


-0.000283468655320 


+2.32371573798 10~ s 


0.942538421086 


5.04507937203 


100 


-4.25340157171-10 17 


-6.46775072494 10~ 49 


0.628752815248 


118.603775377 


500 


-1.16550527223-10 204 


-9.1675098540M0- 240 


2.66350209695 


991.900092992 


1000 


-1.57095384420-10 486 


-2.52129710770-10~ 478 


1.75626461597 


2326.05316169 


2000 


+2.68042467892 -10 1109 


-1.90708173159-10 -955 


10.7685011806 


5351.75953838 


3000 






-2.09002802367 


8617.21920730 



4. Applications and conclusions 
The recurrence formulae Q3.8JI together with (|3.1U|) and (|3,12|) allow in 

principle to compute both rj n and A n with arbitrary accuracy for any value 
of n, but it is clear that with increasing n the number of terms increases 
very rapidly 2 . It would be desirable to simplify the polynomials in and 
(J3.13|) . or at least to reveal some hidden regularities in them, but I doubt 
whether this is possible. The table below demonstrates that it would be 

even impractical to write down explicit expressions for, say, A n for n greater 
than 15 or 20. 



Using the On-Line Encyclopedia of Integer Sequences 
(http://www.research.att.eom/~njas/sequences/l one can see that number of terms 
in rj n is related to the number of partitions of n (partition number, PartitionsP[n] 
in Mathematica notation), which grows like exp (const \/n), whereas the number of 

terms in A n is equal to the number of sums S of positive integers satisfying S < n 
(Sum[PartitionsP[k],{k,l,n}] in Mathematica notation). 
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n 


# of terms m rj n (eq. 13.111) 


# oi terms m \ n (eq. 13.1311 





1 




1 


2 


1 


Z 


3 


3 


3 


5 


6 


4 


7 


11 


5 


11 


18 


10 


56 


138 


20 


792 


2713 


30 


6842 


28628 


40 


44583 


215307 


50 


239943 


1295970 



Using the above formulae (|3.8|h (|3.1()j) and (|3.13|) I have computed quite 
a lot of initial values of rj n and A n . First it was necessary to tabulate 
Stieltjes constants j n with sufficient number of significant digits. In order 
to obtain these I used Mathematica 5 which can handle arbitrary precision 
numbers and performs automatically full control of accuracy in numerical 
calculations. (For details concerning Mathematica interval arithmetic see 
e.g. |12j). This part of computations took over 60 hours on AMD 1667 MHz 
processor. 

Recently Kreminski .8.. published an effective method of computing Stielt- 
jes gamma using Newton-Cotes integration algorithm. His method would 
be of considerable interest since, due to some "hardcoded" limitations, the 
current Mathematica version can't give j n (with sufficient accuracy) beyond 
n « 2050. 

The main calculations (r] n and A n ) were also time consuming (about 20 
hours) and required also considerable amount of computer memory (3x256 
Mb). In particular, having 2000 pre-computed Stieltjes constants, with 800 
significant digits each, I calculated 2000 r\ n and almost 3300 A n . Due to finite 
accuracy and the obvious phenomenon of error accumulation, the number 
of significant digits in r\ n and A n decreases with increasing n (see Fig. 6). 

In order to verify the computations as well as to compare various pack- 
ages I also tried to repeat the whole procedure using Maple 8. However, it 
turned out that it is impossible to get Stieltjes constants j n for relatively 
small n =100 with the required precision 800 significant digits. 

(Place Figure 6 about here.) 

The main conclusion which stems from the above calculations is contained 
in the following plots showing the trend of A (Fig. 7a) and the oscillating 
part of A (Fig. 7b). Their sum gives the coefficients which appear in Li's 
criterion for RH. Note that the scales on both plots differ by nearly two 
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orders of magnitude. As mentioned before, it is easy to show that the trend 
(12.5(1 is strictly growing. Therefore, if the oscillations were bounded or, at 
least, if their amplitude would grow with n slower than the trend, then 
RH would be true. In other words, we have a new RH criterion, which is 
simply a reformulation the original Li's result, but from the viewpoint of 
the present paper it has an obvious interpretation. It states that if for all 
positive integer n 

then RH is true. The numerical data gathered so far and presented in 
Figures 7a and 7b is in its favor. Of course, one should bear in mind that in 
number theory the numerical evidence, no matter how "convincing", may 
be just illusory. In fact, Oesterle observed that if the first n zeta zeros are 
on the critical line, then the Li positivity should hold for about the first 
n 2 Li coefficients (see pQ, p. 441). Therefore, direct numerical search for a 
possible counterexample to RH using Li's criterion is rather a hopeless task. 

Finally I would like to stress out that so far there are no published ex- 
tensive tables of Li's coefficients. Several numerical values of A ra are given 
in . All the numerical data obtained during preparation of this paper as 
well as appropriate Mathematica notebook are available from the author. 

(Place Figures 7a and 7b about here.) 

Acknowledgments. I would like to express my gratitude to Prof. Jef- 
frey C. Lagarias, AT&T Labs, for confirming the effect of tiny oscillations, 
as well as for several remarks concerning further development of the idea 
presented in this paper. I would also like to thank Dr. Mark W. Coffey, 
Department of Physics, Colorado School of Mines, for directing my atten- 
tion to some misprints as well as for sending me a preliminary version of his 
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Note added in proof. After completing the calculations I become aware 
of the paper by Keiper 0. His Figure 1 looks similar to Figure 7b of the 
present paper. However it is not exactly the same, insofar as he is subtract- 
ing off an approximation to the " main term" . But presumably the error is a 
constant plus (9(l/n) term, in view of my asymptotic series approximation 
(|2.8[) . so it must be pretty close. 
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Figure captions 

• Figure 1. Distribution of zeroes of £ in the complex plane. 

• Figure 2. Understanding the Li's lambda: inevitable difficulties en- 
countered when calculating numerically both parts of A n . The se- 
quence of consecutive derivatives is to be taken near s = 0. 

• Figure 3. Mobius transformation of the complex plane used by Li. 
The left fragment of the picture is just Fig. 1 reduced to its essential 
part: the critical strip. The half-plane 5fts > \ (left picture) is 
mapped into the unit disk \z\ < 1 (right picture). 

• Figure 4. Plot of l/|C(xri)l on a sma U part of the transformed 
complex plane containing all nontrivial zeroes. On the right fragment 
of Fig. 2 this part of complex plane is a small, very narrow rectangle 
near z = 1. Nontrivial zeroes are visible as sharp "pins". White 
dots are added to help visualize that the zeroes indeed lie on a circle 
(which looks rather like an ellipse here since the scales on tflz and 
Qz are different). The apparent lack of peaks in the center is an 
artifact. All complex zeroes are very crowded near z = 1 and the 
corresponding peaks are increasingly thinner. Obtaining a better 
picture would require much higher density of points in which values 
of zeta are calculated (hence more computer memory) and much 
higher resolution of the picture. 
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• Figure 5. Signs of the coefficients of matrix c (|3.8j) for k = 100 with 
rows and columns labelled as in (|3.9ft . Little white squares denote 
plus sign, black squares denote minus sign; grey squares mark unused 
entries of the matrix. 

• Figure 6. Accuracies of various numbers used in this paper: 7 n , 

r\ n and A n . Having 2000 precomputed Stieltjes constants j n , with 
800 significant digits each, I could (using Mathematica 5) obtain 

2000 coefficients r\ n and about 3300 oscillating parts of lambda, A n , 
both with linearly decreasing accuracy. The accuracies of r] n and 

A n decrease with n due to complicated error cumulating but almost 
perfect linearity of their dependence is a priori not so obvious because 
the number of terms in ()3.11|) and (|3.13|) grows fast with increasing 
n. In particular, the fact that the accuracy of r] n decreases faster 

than accuracy of A n is rather counter intuitive. 

• Figures 7a and 7b. The trend of A n (7a) in comparison with the 
oscillating part of A ra (7b). Note different vertical scales. In fact, the 
sum of the trend and the oscillating part, i.e. full A n , would look 
exactly like the upper plot since the amplitude of the oscillations is 
smaller than the thickness of the graph line. 

Astronomical Observatory of the Jagiellonian University, Orla 171, 30-244 
Cracow, Poland 
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nontrivial 
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40i 
30i 
20i 
1 0i 



trivial zeros 



14 -12 -10 -8 -6 -4 -2 

-10i 

-20i 

-30i 

-40i 



critical line 



no zeros 



1/2 + i 49.773332... 
1/2 + i 43.005151... 

1/2 + i 43.327073... 
1/2 + i 40.918719... 

1/2 + i 37.586178... 

1/2 + i 32.935062... 
1/2 + i 30.424376... 

1,2 + i 25.010358... 
1/2 + i 21.022040... 

12 + i 14.134725... 

simple pole 

i i_ 

2 4 6 



1/2 - i 14.134725. 

1/2 - i 21.022040. 
1/2 - i 25.010353. 

12 - i 30.424376. 
12 - i 32.935062. 

1/2 - i 37.536173. 

i.2 - i 40.913719. 
1/2 - i 43.327073. 

1/2 - i 43.005151. 
1/2 - i 49.773332. 

critical strip 



° 2i 
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Figure 3 




Figure 4. 
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Figure 7b 



